January,  1995 


LIDS-P  2287 


Research  Supported  By: 

ONRN00014-91-J-1004 
ARPA  F49620-93-1-0604 
AFOSR  F49620-95-1-0083 


Fractal  Estimation  using  Models  on  Multiscale  Trees 


Fieguth,  P.W. 
Willsky,  A.S. 


Report  Documentation  Page 

Form  Approved 

0MB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 

1.  REPORT  DATE 

JAN  1995 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-01-1995  to  00-01-1995 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Fractal  Estimation  using  Models  on  Multiscale  Trees 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Massachusetts  Institute  of  Technology, Laboratory  for  Information  and 
Decision  Systems, 77  Massachusetts  Avenue, Cambridge, MA, 02139 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 

OF  PAGES 

13 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98} 

Prescribed  by  ANSI  Std  Z39-18 


January  1995 


T.IDS-P-2287 


Fractal  Estimation  using  Models  on  Multiscale  Trees 

Paul  W.  Fieguth  Alan  S.  Willsky 


Abstract 

In  this  paper  we  estimate  the  fractal  dimension  of  stochastic  processes  with  l//-like  spectra 
by  applying  a  recently-introduced  multiresolution  framework.  This  framework  admits  an  effi¬ 
cient  likelihood  function  evaluation,  allowing  us  to  compute  the  maximum  likelihood  estimate 
of  this  fractal  parameter  with  relative  ease.  In  addition  to  yielding  results  that  compare  well  to 
other  proposed  methods,  and  in  contrast  to  other  approaches,  our  method  is  directly  applicable, 
with  at  most  very  simple  modification,  in  a  variety  of  other  contexts  including  fractal  estimation 
given  irregularly  sampled  data  or  nonstationary  measurement  noise  and  the  estimation  of  fractal 
parameters  for  2-D  random  fields. 


1  Introduction 


Many  natural  and  human  phenomena  have  been  found  to  possess  l//-like  spectral  properties,  which 
has  led  to  considerable  study  of  1  //  processes.  One  class  of  such  processes  that  is  frequently  used 
because  of  its  analytical  convenience  and  tractability  is  the  class  of  fractional  Brownian  motion 
(fBm)  processes,  introduced  by  Mandelbrot  and  Van  Ness[8].  For  practical  computation  purposes 
we  consider  only  sampled  versions  of  continuous  time  fBm  processes  B{t),  i.e., 

B[k]  =  B{kAt)  k€Z  (1) 
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for  which  the  associated  nonstationary  covariance  is 

E  [B[k],B[m]]  =  -\k- 

where  cr  and  H  are  scalar  parameters  which  completely  characterize  the  process.  H,  the  quan¬ 
tity  which  we  wish  to  estimate,  determines  the  fractal  dimension  D  =  2  —  H  oi  the  process. 
Previous  estimators  have  been  developed  addressing  this  problem,  notably  those  of  Wornell  and 
Oppenheim[ll],  Kaplan  and  Kuo[4],  Tewfik  and  Deriche[lO],  and  Flandrin[3].  The  exact  maximum 
likelihood  (ML)  calculation  for  the  fractal  dimension  of  fBm  is  computationally  difficult  (see  [10]); 
to  address  this  difficulty,  fractal  estimators  typically  fall  into  one  of  the  two  following  classes  to 
achieve  computational  efficiency: 

1.  optimal  algorithms,  admitting  an  efficient  solution,  based  on  a  l//-model  other  than  fBm; 

2.  approximate  or  suboptimal  algorithms  developed  directly  from  the  fBm  model. 

Our  approach  and  that  of  [11]  fall  into  the  former  category,  while  the  methods  in  [3,  4,  9]  fall  into 
the  latter.  In  particular  the  approach  in  [11]  is  based  on  a  1//  process  constructed  using  a  wavelet 
basis  in  which  the  wavelet  coefficients  are  independent,  with  variances  that  vary  geometrically 
with  scale  with  exponent  equal  to  H.  The  method  in  [4]  determines  the  exact  statistics  of  the  Haar 
wavelet  coefficients  of  the  discrete  fractional  Gaussian  noise  (DFGN)  process  F[k]  =  B[k  -f  1]  —  B[k] 
and  then  develops  an  estimator  by  assuming,  with  some  approximation,  that  the  coefficients  are 
uncorrelated. 

The  goal  of  our  research,  on  the  other  hand,  is  the  development  of  a  fast  estimator  for  H  that 
functions  under  a  broader  variety  of  measurement  circumstances,  for  example  the  presence  of  gaps  in 
the  measured  sequence,  measurement  noise  having  a  time-varying  variance  and  higher  dimensional 
processes  (e.g.,  2-D  random  fields).  The  basis  for  accomplishing  this  is  the  utilization  of  a  recently- 
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introduced  multiscale  estimation  framework.  [1.  5]  The  next  section  gives  a  brief  description  of  this 
framework,  followed  by  the  development  of  the  estimator,  and  finally  a  description  of  estimation 
results. 

2  Multiscale  Framework 

In  the  framework  developed  in  [1,  5]  stochastic  models  are  constructed  recursively  in  scale  or  in 
resolution  on  multilevel  trees.  Specifically,  let  s  index  the  nodes  of  a  tree  T  (refer  to  Figure  1) 
which,  for  the  purposes  of  this  paper,  may  be  considered  a  dyadic  tree  although  the  framework 
permits  much  greater  flexibility.  Let  So  G  T  designate  the  root  node  of  T ;  also  let  sj  denote  the 
parent  node  oi  s  ^  So-  Each  node  s  €  T  has  associated  with  it  a  state  vector  2;(s)  and  possibly  an 
observation  vector  y{s).  Stochastic  models  on  T  are  written  recursively  in  this  form 

a;(s)  =  A(s)x(s7)  +  G(s)u;(s)  'is£T,s^So  (3) 

and  where  w{s)  is  a  white  Gaussian  noise  process  with  identity  covariance.  Similarly,  noisy  obser¬ 
vations  of  the  process  are  permitted  on  an  arbitrary  subset  of  the  tree  nodes: 

y{s)  =  C(s)x(s)  +  t)(s)  Vs  €  C?  C  T  (4) 

where  x(s)  is  white  and  Gaussian  with  covariance  R{s).  In  general,  tree  variables  at  all  scales  may 
be  physically  meaningful  and  measurable,  or  they  may  be  abstract  -  a  by-product  of  achieving  the 
desired  statistics  on  the  finest  level  of  the  tree.  In  this  paper  coarse  scale  nodes  are  abstract,  with 
a  l//-like  process  residing  on  the  finest  scale. 

For  those  multiscale  stochastic  models  which  can  be  written  in  the  form  (3), (4)  the  following 
two  problems  possess  extremely  efficient  algorithmic  realizations[l,  6]: 

1.  Given  observations  y{),  determine  the  optimum  least-square  estimate  for  x() 
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2.  Determine  the  likelihood  I 


AQ, GO, CO, ROM) 


of  a  set  of  observations  y() 


The  latter  algorithm  permits  the  estimation  of  any  parameter  embedded  in  the  multiscale  model  by 
maximizing  the  likelihood  function;  this  is  the  very  problem  which  we  have  set  out  to  solve  in  this 
paper:  by  formulating  an  appropriate  multiscale  model  A{s,H),G{s,H),C{s),R{s)  an  estimator 
for  H  may  be  written  abstractly  as 


H  =  arg^  max  I 


A{s,H),G{s,H),C{s),R{s),B[k] 


(5) 


As  was  the  case  with  Wornell  and  Oppenheim[ll],  we  do  not  construct  an  exact  model  of  fBm,  but 
rather  choose  an  appropriate  approximation  -  in  our  case  within  this  multiscale  framework.  The 
selection  of  such  a  multiscale  model  is  achieved  in  the  next  section. 


3  Fractal  Estimator 


The  statistical  self-similarity  of  fBm  makes  the  application  of  wavelets  a  logical  choice.  Kaplan  and 
Kuo[4]  apply  the  Haar  wavelet  to  the  incremental  process  F[fc],  and  Wornell  and  Oppenheim[ll] 
apply  higher  order  Daubechies  wavelets  to  B[fc].  We  use  the  multiscale  framework  just  described 
to  develop  a  Haar  wavelet  multiscale  stochastic  model  which  applies  directly  to  B[k].  This  choice 
of  wavelet  is  motivated  by  the  particularly  simple  realization  of  the  Haar  wavelet  in  our  multiscale 


framework  by  using  a  dyadic  tree  structure  (see  Figure  1): 


Coarse  Scales  ; 


Ix{s)  = 
xis)  = 


1  +1 
0  0 
1  -1 
0  0 


a:(s7)  + 
a:(s7)  + 


g{s,H)w{s) 

g(s,H)w(s) 


s  is  left  descendant 
of  its  parent 
s  is  right  descendant 
of  its  parent 


Finest  Scale  : 


x{s)  =  [  1  +1  j  x(s7)  -k  0  •  w{s) 

x{s)  =  [  1  — 1  ]  a:(s7)  +  0  •  w{s) 
.  y(s)  =  a:(s)  -ku(s) 


s  is  left  descendant 
of  its  parent 
s  is  right  descendant 
of  its  parent 


(6) 


(7) 
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That  is,  at  coarse  scales  x{s)  consists  of  two  scalars:  a  coarse  approximation  to  the  1//  process,  and 
a  detail  coefficient.  The  detail  coefficient  equals  the  difference  in  the  coarse  l//-like  representation 
between  node  s  and  its  two  children,  where  the  sign  of  this  difference  depends  on  the  parity  of 
the  child  (i.e.,  left  vs.  right).  At  the  finest  scale  x(s)  is  a  single  scalar,  representing  a  sample  of  a 
l//-like  process,  and  measurements  of  the  actual  fBm  sequence  appear  as  observations  y{s)  at  the 
finest  scale.  Note  that  the  1/f  process  described  by  (6), (7)  does  not  yield  a  finest  scale  process 
that  is  exactly  an  fBm  process  (and  thus,  as  with  the  technique  in  [11],  our  model  does  not  exactly 
match  the  statistics  of  the  process  to  be  estimated).  However  by  an  appropriate  choice  of  the 
remaining  model  parameters  we  can  produce  a  process  with  the  same  type  of  1//  behavior. 

The  elements  which  remain  to  be  determined  in  the  above  multiscale  model  are  the  g{s,H): 
the  variance  of  the  detail  wavelet  coefficients  between  node  s  and  its  children.  Expressions  for  the 
statistics  of  the  wavelet  decomposition  of  fBm  have  been  determined  by  others  [3,  9],  however  the 
self-statistics  for  the  special  case  of  the  Haar  wavelet  are  easily  computed  as  follows: 

•  Let  Bo[k]  =  B[k]  which  is  the  fBm  process  of  interest. 

•  Define  Bm[k]  as  the  process  obtained  by  coarsening  B[k]  m  times, 

Bm[k]  =  {Bm-l[2k]  +  Bm-l[2k  +  l])/2  (8) 


a  relation  which  follows  from  the  multiscale  model  of  (6). 


From  (8)  it  follows  that 


Bm-i[k] 


B  [2'^~^k  +  i] 

i=0 


-1,  .1  i  -  1 


Bm^,[k  +  l\~Bm^,[k]  =  E  + 

i=0 

2(2'^-! -1) 

=  E  F[2'^-^k  +  i\ci 


i=0 


(9) 

(10) 

(11) 
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(12) 


•  From  the  stationarity  of  the  increments  process  F[fc],  and  from  the  symmetry 

Bm[k]  T  Bm-i[2k]  =  -  {Bm[k]  -  B,n-i[2k  +  1]} 

(11)  reduces  to  the  desired  variance  expression 

2(2”*-'-!)  2(2'"-l-l)+min(0,-i) 

E[{B^[k]-B^_i[2k])^]  =  -  ^fH]  E  CiCi+j=^g\s,H){U) 

t=— 2(2'"“! -1)  j=-min(0,i) 

where  s  is  any  node  on  scale  m  of  the  tree,  and  where  Ap  is  the  covariance  function  of  F[k]: 

\w  =  y[ii+ir+ii-ii“-2iir]  (14) 

By  way  of  comparison,  it  has  been  shown[ll]  that  1//  processes,  of  which  fBm  is  a  subset,  may  be 

approximated  by  wavelet  synthesis  in  which  the  wavelet  coefficients  vary  exponentially  with  scale: 

gl  =  i.e.,  log2  ^  =  i?  (15) 

9m— l 

Table  1  shows  the  scale  to  scale  variance  ratios  as  predicted  by  (13).  The  deviation  from  the 
approximate  scaling  law  of  (15)  is  most  pronounced  at  low  iJ;  it  is  this  deviation  which  leads  to  a 
bias  for  those  estimators  based  on  (15),  as  shown  in  Table  2. 

The  actual  estimator  for  H,  based  on  the  multiscale  model  of  (6), (7),  takes  precisely  the  form 
as  outlined  in  (5),  in  which  the  likelihood  maximization  is  performed  using  standard  nonlinear 
techniques  (e.g.,  the  section  search  method  of  Matlab). 

4  Experimental  Results 

Sixty  four  fBm  sample  paths,  each  having  a  length  of  2048  samples,  were  generated  using  the 
Cholesky  decomposition  method  of  [7],  precisely  the  same  approach  as  in  Kaplan  and  Kuo[4]  whose 
experimental  results  form  the  basis  of  comparison  with  ours. 
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The  performance  of  three  fBm  estimators  is  compared  in  Table  3.  The  bias  in  the  estimator  of 
[11]  for  low  H,  as  was  argued  earlier  based  on  Table  1,  is  evident.  Also  recall  that  the  multiscale 
model  of  (6), (7)  assumed  the  wavelet  detail  coefficients  to  be  uncorrelated;  this  assumption  becomes 
progressively  poorer  as  H  increases [3],  leading  to  an  increase  in  the  error  variance  for  our  estimator 
at  large  H.  Nevertheless,  our  method  still  performs  reasonably  well  over  quite  a  wide  range  of  values 
of  H.  Moreover,  using  the  techniques  developed  in  [5],  we  can  construct  higher-order  multiscale 
models  which  account  for  most  of  the  residual  correlation  in  the  wavelet  coefficients.  The  same 
likelihood  procedure  applied  to  these  higher-order  models  would  then  yield  even  better  results, 
closely  approaching  the  exact  ML  estimator  based  on  the  exact  fBm  statistics.  However  since  fBm 
itself  is  an  idealization,  the  benefit  in  practice  of  such  higher-order  models  over  that  based  on  the 
low-order  model  (6), (7)  depends  upon  the  application. 

Finally,  as  we  have  said,  our  approach  applies  equally  well  in  a  variety  of  other  settings.  For 
example.  Table  4  illustrates  the  performance  of  our  estimator  under  non-uniform  sampling  (by 
randomly  discarding  10%  of  the  measurements),  and  under  a  varying  measurement  noise  variance. 
Both  of  these  special  cases  are  accomplished  with  essentially  no  change  in  the  algorithm.  In 
addition,  by  using  a  quadtree  rather  than  a  dyadic  tree  we  can  also  apply  these  techniques  in  2-D. 
An  example  of  such  an  application  to  the  estimation  of  non-isotropic  fractal  parameters  for  a  2-D 
random  field  based  on  irregular,  nonstationary  data  is  given  in  [2], 
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List  of  Captions; 


Table  1: 


Table  2; 


Table  3: 


Table  4: 


Fig.  1: 


Haar  wavelet  coefficient  standard-deviation  ratios  as  a  function  of  H  and  scale:  represents 

the  wavelet  coefficient  variance  at  scale  i,  where  the  finest  scale  is  i  =  0.  The  deviation  of  the 
variance  progression  from  an  exponential  law  is  most  pronounced  at  fine  scales  and  for  low 
values  of  H. 

The  estimator  in  the  top  row  is  based  on  the  premise  that  wavelet  coefficient  variances  are 
exponentially  distributed,  leading  to  a  biased  estimate  H. 

Estimation  results  for  three  estimators,  based  on  64  fBm  sample  paths,  each  of  length  2048 
samples,  with  no  additive  noise.  The  experimental  results  for  the  first  two  estimators  are 
from  [4], 

Performance  of  the  fBm  estimator  for  two  examples:  irregular  sampling  (removal,  at  random, 
of  10%  of  the  measurements),  and  nonstationary  measurement  noise  variance  (noise  standard 
deviation  cr[^]  =  ^exp{  — [(A:  —  1024)/500]^}).  In  both  cases  the  results  are  based  on  64  fBm 
sample  paths,  each  of  length  2048  samples. 

Dyadic  tree  structure  used  for  the  estimator  of  this  paper. 


9 


Variance  Ratio: 

77  =  0.25  iJ--0.50  iJ  =  0.75  H  =  0.9 

iog2 

log2  (pr/pe) 
log2  (^e/ps) 
log2  (55/54) 
log2  (54/53) 
log2 (53/52) 
log2  (52/51) 
log2  (51/50) 

0.250 

0.249  i 

0.247  0.500  : 

0.242  0.499  0.750  i 

0.228  0.496  0.749  0.900 

0.188  0.484  0.745  0.898 

0.091  0.437  0.727  0.892 

-0.084  0.292  0.650  0.861 

Table  1: 


Variance  Rule: 

H  =  0.25 

H  =  0.50 

H  =  0.75 

H  =  0.90 

Variances  assumed  ex¬ 
ponential  with  scale 

H: 

0.05 

0.40 

0.70 

0.91 

Variances  based  on  ex¬ 
act  result 

H: 

0.24 

0.51 

0.75 

0.92 

Table  2: 
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Estimator 

H  =  0.25 

H  =  0.50 

H  =  0.75 

H  =  0.90 

H 

0.082 

0.398 

0.683 

0.846 

W.O. 

a. 

H 

0.022 

0.021 

0.026 

0.021 

(H  -  HU, 

0.169 

0.109 

0.072 

0.058 

H 

0.252 

0.499 

0.748 

0.899 

K.K. 

0.017 

0.017 

0.017 

0.017 

(H  -  HUs 

0.017 

0.017 

0.017 

0.017 

Multiscale 

H 

0.249 

0.503 

0.768 

0.919 

Haar 

0.011 

0.019 

0.050 

0.109 

{H  -  HUs 

0.011 

0.019 

0.054 

0.110 

Table  3: 


Circumstance 

H  =  0.25 

H  =  0.50 

H  =  0.75 

H  =  0.90 

Irregular 

H 

0.246 

0.507 

0.781 

0.937 

Sampling 

0.033 

0.044 

0.076 

0.124 

(H  -  HUs 

0.033 

0.045 

0.082 

0.128 

Nonstationary 

H 

0.268 

0.511 

0.769 

0.918 

Measurement 

0.011 

0.019 

0.051 

0.109 

Noise  Variance 

{H  -  HUs 

0.021 

0.022 

0.054 

0.110 

Table  4: 
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